\documentclass{article}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsthm}
\usepackage{mathtools}
\usepackage{hyperref}
\DeclarePairedDelimiter\abs{\lvert}{\rvert}
\newtheorem{proposition}{Proposition}
\def\E{\mathbb{E}}
\def\R{\mathbb{R}}
\DeclareMathOperator*{\tr}{tr}
\title{High order moment of uniform distribution on the sphere}
\author{Feng Zhao}
\begin{document}
\maketitle
Consider $X_1, \dots, X_n \sim N(0, 1)$.
Let $X = \sqrt{X_1^2 + \dots + X_n^2}$.
$Y_i = \frac{X_i}{X}$.
Then $Y=(Y_1, \dots, Y_n)$ are uniformly distributed on the sphere.
Using symmetric property we have $\E[Y_i^2] = \frac{1}{n}$.
Our problem is to compute the higher order moments of $\E[Y_i^k]$
where $k$ is even. Notice that $\E[Y_i^k]=0$ if $k$ is odd.

First we consider the simple case in which $n=2$ and we show that
\begin{equation}\label{eq:n2}
\int_{\R}\int_{\R} \frac{x^{2k}}{(x^2+y^2)^k}
    p(x)p(y)dxdy = \frac{(2k-1)!!}{(2k)!!}
\end{equation}
where $p(x) = \frac{1}{\sqrt{2\pi}} \exp(-\frac{x^2}{2})$.
We can let $y=xt$ and the integral is transformed to
\begin{align*}
\int_{\R}\int_{\R} \frac{x^{2k}}{(x^2+y^2)^k} p(x)p(y)dxdy & =
\frac{2}{\pi}\int_{0}^{+\infty}\int_{0}^{+\infty} \frac{x}{(1+t^2)^k}
\exp(-\frac{1}{2}x^2(1+t^2))dxdt \\
&=
\frac{2}{\pi}\int_{0}^{+\infty} \frac{1}{(1+t^2)^{k+1}}dt \\
&= \frac{2}{\pi}\int_{0}^{\frac{\pi}{2}} cos^{2k}(u)dt \\
\end{align*}
The integral of $sin(x)$ over $[0, \frac{\pi}{2}]$ is a well-known
\href{https://math.stackexchange.com/questions/50447/
integration-of-powers-of-the-sin-x}{formula}.
And we conclude that the right hand size of equation \eqref{eq:n2} equals
$\frac{(2k-1)!!}{(2k)!!}$.

Next we consider the case when $n>2$.
Let $Z = X_2^2 + \dots + X_{n-1}^2$ and $Z$ is chi-square distribution.
We have
\begin{equation}\label{eq:n3}
\int_{\R}\int_{0}^{+\infty} \frac{x^{2k}}{(x^2+z)^k} p(x)q(z)dxdz
= \frac{\Gamma(\frac{n}{2}) \Gamma(\frac{2k+1}{2})}
{\sqrt{\pi} \Gamma(\frac{2k+n}{2})} = \prod_{i=0}^{k-1} \frac{1+2i}{n+2i}
\end{equation}
where $q(z) = \frac{1}{2^{\frac{n-1}{2}}\Gamma(\frac{n-1}{2})}
z ^{\frac{n-1}{2} - 1} \exp(-\frac{z}{2})$
This time we let $z=x^2 t$ and the integral is transformed to
\begin{align*}
\int_{\R}\int_{0}^{+\infty} \frac{x^{2k}}{(x^2+z)^k} p(x)q(z)dxdz & = \\
\frac{1}{2^{\frac{n}{2}-1}\sqrt{\pi}\Gamma(\frac{n-1}{2})}
\int_{0}^{+\infty}t^{\frac{n-1}{2}-1}
(1+t)^{-k}dt\int_{0}^{+\infty}x^{n-1} \exp(-\frac{1}{2}x^2(1+t)) dx &=
\\
\frac{1}{2^{\frac{n}{2}-1}\sqrt{\pi}\Gamma(\frac{n-1}{2})}
\int_{0}^{+\infty}t^{\frac{n-1}{2}-1}(1+t)^{-k-\frac{n}{2}}dt
\int_{0}^{+\infty}x^{n-1} \exp(-\frac{1}{2}x^2) dx &=
\\
\frac{2}{\sqrt{\pi}\Gamma(\frac{n-1}{2})}
\int_{0}^{+\infty}t^{\frac{n-1}{2}-1}(1+t)^{-k-\frac{n}{2}}dt
\int_{0}^{+\infty}x^{n-1} \exp(-x^2) dx &=
\\
\frac{\Gamma(\frac{n}{2})}{\sqrt{\pi}\Gamma(\frac{n-1}{2})}
\int_{0}^{+\infty}t^{\frac{n-1}{2}-1}(1+t)^{-k-\frac{n}{2}}dt
&= \\ (t=\tan^2 x)
\frac{\Gamma(\frac{n}{2})}{\sqrt{\pi}\Gamma(\frac{n-1}{2})}
2 \int_0^{\frac{\pi}{2}} (\sin x)^{n-2} (\cos x)^{2k} dx
&= \\ \frac{\Gamma(\frac{n}{2})}{\sqrt{\pi}\Gamma(\frac{n-1}{2})}
B(\frac{n-1}{2}, \frac{2k+1}{2})
&= \\ \frac{\Gamma(\frac{n}{2})\Gamma(\frac{2k+1}{2})}
{\sqrt{\pi} \Gamma(\frac{2k+n}{2})}
\end{align*}

Actually, we can caculate the exact distribution of $Y_1$.
Suppose $w(n, q) = \frac{1}{2^{nq/2} \Gamma_q(\frac{n}{2})}$ which is called Wishart constant. $\Gamma_q$ is the multivariate Gamma function.
The Proposition 7.3, Chapter 7 of \cite{eaton1989group} says:
\begin{proposition}\label{prop:dens}
Suppose $X$ is a random orthogonal matrix, $Y$ is the $p\times q$ upper submatrix of $X$. $q\leq p $ and $p+q \leq n$. Then the density function of $Y$ is given by

\begin{equation}
f(y) = \frac{w(n-p, q)}{\sqrt{2\pi}^{pq}w(n,q)} \abs{I_q - y^T y}^{\frac{n-p-q-1}{2}}I_0(y^T y)
\end{equation}
where $I_0$ is an indicator function of the $q\times q $ symmetric matrices. $I_0(A)=1$ if all eigenvalues of $A$ falls within $(0,1)$. Otherwise $I_0(A)=0$.
\end{proposition}
Using this proposition we can get the density function of $Y_1$ as
\begin{equation}
p(y) = \frac{\Gamma(\frac{n}{2})}{\sqrt{\pi}\Gamma(\frac{n-1}{2})} (1-y^2)^{\frac{n-3}{2}}
\end{equation}
Using the property of beta function 
$\int_{-1}^{1} y^{2k} p(y)dy = \frac{\Gamma(\frac{n}{2})\Gamma(\frac{2k+1}{2})}
{\sqrt{\pi} \Gamma(\frac{2k+n}{2})}$, which is same with Equation \eqref{eq:n3}.
Notice that $\hat{Y}_1 = Y_1^2 (-1 \leq Y_1 \leq 1)$ is beta distribution $B(\frac{1}{2}, \frac{n-1}{2})$. The high order moment of $\hat{Y}_1$ is a well known result.
See \href{https://en.wikipedia.org/wiki/Beta_distribution#Higher_moments}{Higher moments of Beta distribution}.
We also consider the case when the moments of
$\hat{Y} = Y_1^2 + \dots + Y_r^2$ is required.
Let $Z = X_{r+1}^2 + \dots + X_n^2$. We have
\begin{equation}\label{eq:n4}
\int_{0}^{+\infty}\int_{0}^{+\infty} \frac{x^t}{(x+z)^t} p(x)q(z) dxdz =
\frac{\Gamma(\frac{2t+r}{2})\Gamma(\frac{n}{2})}
{\Gamma(\frac{2t+n}{2})\Gamma(\frac{r}{2})} =
\prod_{i=0}^{t-1} \frac{2i+r}{2i+n}
\end{equation}
where $p(x)$ is pdf of chi square distribution with degree $r$ and $q(z)$
is pdf of chi square distribution with degree $n-r$.
\begin{align*}
\int_{0}^{+\infty}\int_{0}^{+\infty} \frac{x^t}{(x+z)^t} p(x)q(z) dxdz =
& \\ \frac{1}{2^{n/2} \Gamma(\frac{n-r}{2})\Gamma(\frac{r}{2})}
\int_{0}^{+\infty}\int_{0}^{+\infty}
\left(\frac{x}{x+z}\right)^t
x^{\frac{r}{2}-1}\exp(-\frac{x}{2})
z^{\frac{n-r}{2}-1}\exp(-\frac{z}{2})dxdz =
& \\ (z=xy) \frac{1}{2^{n/2}\Gamma(\frac{n-r}{2})\Gamma(\frac{r}{2})}
\int_{0}^{+\infty}\int_{0}^{+\infty}
\left(\frac{1}{1+y}\right)^t
x^{\frac{n}{2}-1}\exp(-\frac{x(1+y)}{2})y^{\frac{n-r}{2}-1}dxdy =
& \\  \frac{1}{2^{n/2} \Gamma(\frac{n-r}{2})\Gamma(\frac{r}{2})}
\int_{0}^{+\infty}\left(\frac{1}{1+y}\right)^t
y^{\frac{n-r}{2}-1} dy
\int_{0}^{+\infty} x^{\frac{n}{2}-1}\exp(-\frac{x(1+y)}{2})dx = & \\
\frac{\Gamma(\frac{n}{2})}{\Gamma(\frac{n-r}{2})\Gamma(\frac{r}{2})}
\int_{0}^{+\infty}\left(\frac{1}{1+y}\right)^{t+\frac{n}{2}}
y^{\frac{n-r}{2}-1} dy =
& \\  \frac{\Gamma(\frac{n}{2})}{\Gamma(\frac{n-r}{2})\Gamma(\frac{r}{2})}
B(\frac{n-r}{2}, \frac{2t+r}{2}) =
& \\ \frac{\Gamma(\frac{2t+r}{2})\Gamma(\frac{n}{2})}
{\Gamma(\frac{2t+n}{2})\Gamma(\frac{r}{2})}
\end{align*}

From Proposition \ref{prop:dens}, let $p=r, q=1$. The joint distribution of $Y_1, \dots, Y_r$ can be written as:
$$
p(y_1, \dots, y_r) = \frac{\Gamma(\frac{n}{2})}{\sqrt{\pi}^r \Gamma(\frac{n-r}{2})}(1- y_1^2 - \dots - y_r^2)^{\frac{n-r}{2} -1}
$$
To compute $\E[(Y_1^2 + \dots Y_r^2)^t]$, we
make the spherical transformation and integrate over all those angle components:
\begin{align*}
\frac{\Gamma(\frac{n}{2})}{\sqrt{\pi}^r\Gamma(\frac{n-r}{2})}  \frac{r\pi^{r/2}}{\Gamma(\frac{r}{2}+1)}\int_{0}^1 y^{2t+r-1}(1- y^2)^{\frac{n-r}{2} -1} dy =
\frac{\Gamma(\frac{2t+r}{2})\Gamma(\frac{n}{2})}
{\Gamma(\frac{2t+n}{2})\Gamma(\frac{r}{2})}
\end{align*}
The result is the same with Equation \eqref{eq:n4}. We can further show that $\hat{Y}$ is beta distribution with $B(\frac{r}{2}, \frac{n-r}{2})$ using the techniques similar to \href{https://en.wikipedia.org/wiki/Proofs_related_to_chi-squared_distribution#Derivation_of_the_pdf_for_k_degrees_of_freedom}{Derivation of pdf for chi-square distributions}. Then the higher order moment of $\hat{Y}$ follows naturally. The distribution of $\hat{Y}$ can also be got from Proposition 7.2 in \cite{eaton1989group} by setting $p=r, q=1$.


We like to compute $\E[X_1^{2k_1} X_2^{2k_2}]$. Let $ Z = X_3^2 + \dots + X_n^2$ which is a chi square distribution with $n-2$ degree of freedom.
Then we have
\begin{align}\label{eq:cor}
&\int_{x,y\in \mathbb{R}}\int_{z\geq 0} \frac{x^{2k_1} y^{2k_2}}{(x^2+y^2+z)^{k_1 + k_2}} p(x)p(y)q(z) dxdydz \notag\\
&=\frac{1}{n-2}\frac{(2k_1 -1)!!(2k_2-1)!! (n-3)!!}{(n-2 + 2(k_1 + k_2))!!}
\end{align}
We let $z=(x^2+y^2)t$ first.
\begin{align*}
\int_{x,y\in \mathbb{R}}\int_{z\geq 0} \frac{x^{2k_1} y^{2k_2}}{(x^2+y^2+z)^{k_1 + k_2}} p(x)p(y)q(z) dxdydz & =\\
\frac{1}{2^{\frac{n-2}{2}}\Gamma(\frac{n-2}{2})}\frac{1}{2\pi} \int_{x,y\in \mathbb{R}}\int_{t\geq 0}\frac{x^{2k_1} y^{2k_2}t^{-2+\frac{n}{2}}}{(x^2+y^2)^{k_1 + k_2+1-\frac{n}{2}}(1+t)^{k_1 + k_2}}\exp(-\frac{1}{2}(x^2+y^2)(1+t))dxdydt =
\end{align*}
Secondly we make the transformation $y^2=x^2u$ with $u>0$.
\begin{align*}
\frac{1}{2^{\frac{n-2}{2}}\Gamma(\frac{n-2}{2})}\frac{2}{2\pi}\int_{x,u,t\geq 0}\frac{x^{n-1} u^{k_2 - \frac{1}{2}}t^{-2+\frac{n}{2}}}{(1+u)^{k_1 + k_2+1-\frac{n}{2}}(1+t)^{k_1 + k_2}}\exp(-\frac{1}{2}x^2(1+u)(1+t))dxdudt =
\end{align*}
We can first integrate about $x$ first (we assume the variance $\sigma^2 = \frac{1}{(1+u)(1+t)}$)
and get
\begin{align*}
\frac{2^{n/2} \Gamma(\frac{n}{2})}{2^{\frac{n-2}{2}}\Gamma(\frac{n-2}{2})}\frac{1}{2\pi}\int_{u,t\geq 0} \frac{u^{k_2 - \frac{1}{2}}t^{-2+\frac{n}{2}}\sigma^n}{(1+u)^{k_1 + k_2+1-\frac{n}{2}}(1+t)^{k_1 + k_2}}dudt = \\
\frac{n-2}{2\pi}\int_{u,t\geq 0} \frac{u^{k_2 - \frac{1}{2}}t^{-2+\frac{n}{2}}}{(1+u)^{k_1 + k_2+1}(1+t)^{k_1 + k_2 + \frac{n}{2}}}dudt = \\
\frac{n-2}{2\pi} B(k_2 + \frac{1}{2}, k_1 + \frac{1}{2}) B(-1+\frac{n}{2}, k_1 + k_2 + 1) =\\
\frac{n-2}{2\pi} \frac{\Gamma(k_1 + \frac{1}{2}) \Gamma(k_2 + \frac{1}{2})\Gamma(\frac{n}{2}-1)}{\Gamma(k_1 + k_2 + \frac{n}{2})}
\end{align*}
From Equation \eqref{eq:cor}, we choose $k_1 = k_2 = k$ and we have
$\E[X_1^{2k}X_2^{2k}]=\frac{1}{n-2}\frac{((2k-1)!!)^2}{\prod_{i=0}^{2k} (n-2 + 2i)}$

We consider the decomposition of a $2\times 2$ symmetric matrix:
$$
\begin{pmatrix}
x & y \\
y & z 
\end{pmatrix} = 
\begin{pmatrix}
\cos \theta & -\sin \theta \\
\sin \theta & \cos \theta 
\end{pmatrix}
\begin{pmatrix}
\lambda_1 & 0 \\
0 & \lambda_2
\end{pmatrix} 
\begin{pmatrix}
\cos \theta & \sin \theta \\
-\sin \theta & \cos \theta 
\end{pmatrix}
$$
We then have
\begin{align*}
x & = \cos^2\theta \lambda_1 + \sin^2 \theta \lambda_2 \\
y & = \cos\theta \sin\theta (\lambda_1 - \lambda_2) \\
z & = \sin^2\theta \lambda_1 + \cos^2 \theta \lambda_2
\end{align*}
If we have an integral over all $2\times 2$ positive definite matrices, we can convert the integral to $\theta, \lambda_1, \lambda_2$. The determinant of the Jacobi matrix is $\abs{\lambda_1 - \lambda_2}$. The domain for $(\theta, \lambda_1, \lambda_2)$ is
$\lambda_1 \geq \lambda_2 \geq 0$ and $\theta \in [0,\pi]$.
We show that $\Gamma_2(a) = \sqrt{\pi}\Gamma(a) \Gamma(a-\frac{1}{2})$ where
$\Gamma_2(a) = \int_{S>0} \exp(-\tr(S)) \abs{S}^{a-\frac{3}{2}} dS$.

First we transform the integral to $(\lambda_1, \lambda_2, \theta)$ space and integrate out $\theta$:
$$
\Gamma_2(a) = \pi \int_{\lambda_1 \geq \lambda_2} (\lambda_1 -\lambda_2) \exp(-\lambda_1 - \lambda_2) (\lambda_1 \lambda_2)^{a-\frac{3}{2}} d\lambda_1 d\lambda_2
$$
Next we make another transformation: $\lambda_1 = r \cos^2 \phi, \lambda_2 = r \sin^2 \phi$. The determinant of Jacobi is $2r\sin\phi \cos\phi$ where $ r >0, \phi \in [0, \frac{\pi}{4}]$.
Then we have
\begin{align*}
\Gamma_2(a) =&  2\pi \int_0^{\frac{\pi}{4}} (\cos^2\phi - \sin^2 \phi) (\cos\phi\sin\phi)^{2a-2} d\phi \int_0^{+\infty} r^{2a-1} \exp(-r) dr \\
= & 2\pi \Gamma(2a) \int_0^{\frac{pi}{4}} \cos(2\phi) \frac{1}{2^{2a-2}}\sin^{2a-2} (2\phi)d\phi \\
= & \pi\frac{\Gamma(2a)}{2^{2a-2}} \int_0^1 x^{2a-2}dx \\
= & \frac{\pi\Gamma(2a-1)}{2^{2a-2}}
\end{align*}
Using Legendre duplication formula $\Gamma(z) \Gamma(z+\frac{1}{2}) = 2^{1-2z} \sqrt{\pi} \Gamma(2z)$, let $ z = a - \frac{1}{2} $ we have
$\Gamma_2(a) = \sqrt{\pi}\Gamma(a)\Gamma(a-\frac{1}{2})$.

For $p \times p$ symmetric matrix, we define the multivariate Beta function
as 
\begin{equation}
B_p(a, b) = \int_{X \in S_{p,p}^{+} }\abs{I-X}^{a-\frac{p+1}{2}} \abs{X}^{b-\frac{p+1}{2}}dX
\end{equation}
The integration is about $p \times p$ positive definite matrix. We will show that
\begin{equation}
B_p(a,b) = \frac{\Gamma_p(a)\Gamma_p(b)}{\Gamma_p(a+b)}
\end{equation}
where $\Gamma_p(a)$ is defined as:
\begin{equation}
\Gamma_p(a) =  \int_{X \in S_{p,p}^{+} }\abs{X}^{a-\frac{p+1}{2}} \exp(-tr(X))dX
\end{equation}
The proof is similar with \href{https://en.wikipedia.org/wiki/Beta_function#Relationship_between_gamma_function_and_beta_function}{one dimensional case}. 
First we consider 
$\Gamma_p(a)\Gamma_p(b) = \int_{X,Y \in S_{p,p}^{+} }\abs{X}^{a-\frac{p+1}{2}} \abs{Y}^{b-\frac{p+1}{2}}\exp(-tr(X+Y))dXdY$
Then we make the transformation $X = L(I-T)L^T, Y=LTL^T$ where $Z=LL^T$. We actually transform  the integral to $(T,Z)$ space where $I-T, T,Z$ are positive definite.
The Jacobi of this transformation is $|Z|^{\frac{p+1}{2}}$ from the Proposition 5.11, \cite{eaton}. Therefore, the integral transforms to
$$
\int_{Z,T, I-T \in S_{p,p}^{+} }\abs{I-T}^{a-\frac{p+1}{2}} \abs{T}^{b-\frac{p+1}{2}}\abs{Z}^{a+b-\frac{p+1}{2}}\exp(-tr(Z))dZdT
$$ which equals $\Gamma_p(a+b)B(a,b)$ by definition.

\bibliographystyle{plain}
\bibliography{exportlist}

\end{document}
